Download US Census data on natural resource dependent jobs (total # and median income) within counties and summarize across NEP boundaries:
nep_sf <- st_read(dsn = "./data-raw", layer = "NEP_Boundaries10162018", quiet = TRUE) %>%
mutate(NEP_NAME = recode(NEP_NAME, "San Franciso Estuary Partnership" = "San Francisco Estuary Partnership"))
nep_centroid <- st_centroid(nep_sf)
states <- c("AL", "AZ", "AR", "CA", "CO", "CT", "DE", "FL", "GA",
"ID", "IL", "IN", "IA", "KS", "KY", "LA", "ME", "MD",
"MA", "MI", "MN", "MS", "MO", "MT", "NE", "NV", "NH", "NJ",
"NM", "NY", "NC", "ND", "OH", "OK", "OR", "PA", "RI", "SC",
"SD", "TN", "TX", "UT", "VT", "VA", "WA", "WV", "WI", "WY",
"PR")
metrics <- load_variables(2017, "acs5", cache = TRUE) %>%
filter(grepl("fishing|hunting", label))
years <- lst(2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017)
totaljobs_sf <- get_acs(geography = "county", variables = c(Ag_For_Fish_Hunt = "S2403_C01_003"), #C24050_002
state = states, year = 2017, output = "tidy",
geometry = TRUE)
nr_medinc_sf <- get_acs(geography = "county", variables = c(Med_Income = "S2413_C01_003"), #B24031_003
state = states, geometry = TRUE)
ustotaljobs_sf <- st_transform(totaljobs_sf, st_crs(nep_sf))
usnrmedinc_sf <- st_transform(nr_medinc_sf, st_crs(nep_sf))
nep_jobs_intersects <- st_intersects(nep_sf, ustotaljobs_sf)
nep_medinc_intersects <- st_intersects(nep_sf, nr_medinc_sf)
nep_sel_sf <- ustotaljobs_sf[unlist(nep_jobs_intersects),]
nep_sel2_sf <- nr_medinc_sf[unlist(nep_medinc_intersects),]
nep_jobs <- st_join(nep_sf, nep_sel_sf, join = st_intersects) %>%
sf::st_buffer(dist = 0)
nep_nrmedinc <- st_join(nep_sf, nep_sel2_sf, join = st_intersects) %>%
sf::st_buffer(dist = 0)
nep_jobs_sum <- nep_jobs %>%
select(NEP_NAME, estimate) %>%
group_by(NEP_NAME) %>%
summarise(jobs = sum(estimate))
nep_medinc_sum <- nep_nrmedinc %>%
select(NEP_NAME, estimate) %>%
group_by(NEP_NAME) %>%
summarise(medinc = median(estimate))
Plot US Census natural resource dependent jobs data within the NEP boundaries:
pal <- colorNumeric(palette = "viridis", domain = nep_jobs_sum$jobs, n = 10)
nep_jobs_sum %>%
st_transform(crs = "+init=epsg:4326") %>%
leaflet(width = "100%") %>%
addProviderTiles(provider = "CartoDB.Positron") %>%
addPolygons(popup = ~ paste(NEP_NAME,"</br>","Jobs = ",jobs),
stroke = FALSE,
smoothFactor = 0,
fillOpacity = 0.7,
color = ~ pal(jobs)) %>%
addLegend("bottomright",
pal = pal,
values = ~ jobs,
title = "Natural Resource Dependent Jobs</br>(Farming,Forestry,Fishing & Hunting)",
opacity = 1)
Underlying natural resource dependent jobs estimated within NEP Project Areas:
nep_jobs2 <- nep_jobs_sum %>%
select(NEP_NAME, jobs) %>%
st_set_geometry(NULL) %>%
rename("NEP" = NEP_NAME, "Total Jobs (2017)" = jobs)
nep_jobs2
Plot US Census natural resource dependent jobs, median income data within the NEP boundaries:
pal2 <- colorNumeric(palette = "viridis", domain = nep_medinc_sum$medinc, n = 10)
nep_medinc_sum %>%
st_transform(crs = "+init=epsg:4326") %>%
leaflet(width = "100%") %>%
addProviderTiles(provider = "CartoDB.Positron") %>%
addPolygons(popup = ~ paste(NEP_NAME,"</br>","Median Income = ",medinc),
stroke = FALSE,
smoothFactor = 0,
fillOpacity = 0.7,
color = ~ pal2(medinc)) %>%
addLegend("bottomright",
pal = pal2,
values = ~ medinc,
title = "Natural Resource Dependent Jobs</br>(Median Annual Income)",
opacity = 1)
Underlying natural resource dependent jobs, median incomes estimated within NEP Project Areas:
nep_medinc2 <- nep_medinc_sum %>%
select(NEP_NAME, medinc) %>%
st_set_geometry(NULL) %>%
rename("NEP" = NEP_NAME, "Median Income ($, 2017)" = medinc)
nep_medinc2
Import US BEA data on natural resource dependent jobs (total county income) and summarize across NEP boundaries:
userSpecList <- list('UserID' = beaKey, #Need a US BEA API key here
'Method' = 'GetData',
'datasetname' = 'Regional',
'TableName' = 'CAINC5N',
'LineCode' = '100',
'GeoFIPS' = 'COUNTY',
'Year' = '2017')
nr_income <- beaGet(userSpecList, asTable = TRUE) %>%
filter(str_detect(GeoFips, "....0")==FALSE) %>%
mutate(GEOID = as.character(GeoFips))
nr_income$GEOID <- gsub(" ", "", nr_income$GEOID, fixed = TRUE)
nr_inc_sf <- left_join(nep_sel_sf, nr_income, by = c('GEOID' = 'GEOID'))
nep_totinc <- st_join(nep_sf, nr_inc_sf, join = st_intersects) %>%
sf::st_buffer(dist = 0)
nep_inc_sum <- nep_totinc %>%
select(NEP_NAME, DataValue_2017) %>%
group_by(NEP_NAME) %>%
summarise(totinc = sum(DataValue_2017, na.rm = TRUE))
nep_inc_sum$totinc <- ifelse(nep_inc_sum$totinc==0,NA,nep_inc_sum$totinc*1000)
Plot US BEA natural resource dependent jobs total income estimates within the NEP boundaries (where available):
pal3 <- colorNumeric(palette = "viridis", domain = nep_inc_sum$totinc, n = 10)
nep_inc_sum %>%
st_transform(crs = "+init=epsg:4326") %>%
leaflet(width = "100%") %>%
addProviderTiles(provider = "CartoDB.Positron") %>%
addPolygons(popup = ~ paste(NEP_NAME,"</br>","Total Income (2017) = ",totinc),
stroke = FALSE,
smoothFactor = 0,
fillOpacity = 0.7,
color = ~ pal3(totinc)) %>%
addLegend("bottomright",
pal = pal3,
values = ~ totinc,
title = "Natural Resource Dependent Jobs</br>(Total Annual Income, if reported)",
opacity = 1)
Underlying natural resource dependent jobs, total personal income estimated within NEP Project Areas:
nep_inc2 <- nep_inc_sum %>%
select(NEP_NAME, totinc) %>%
st_set_geometry(NULL) %>%
rename("NEP" = NEP_NAME, "Total Personal Income ($, 2017)" = totinc)
nep_inc2
Import GPRA data related to habitat restoration activities across the NEP Project areas, and produce animated map of restoration activities by year. These restoration activities support the natural resource dependent jobs described above.
habitat <- read.csv('./data-raw/GPRA_NEP_HabCat_by_activity_.csv', header=TRUE)
hab_sum <- habitat %>%
mutate(Activity = recode(Activity, "-" = "Unspecified")) %>%
mutate(NEP.Name = recode(NEP.Name, "Santa Monica Bay Restoration Project" = "Santa Monica Bay Restoration Foundation",
"San Francisco Estuary Project" = "San Francisco Estuary Partnership",
"Lower Columbia River Estuary" = "Lower Columbia Estuary Partnership",
"Mobile Bay Estuary Program" = "Mobile Bay National Estuary Program",
"Indian River Lagoon NEP" = "Indian River Lagoon National Estuary Program",
"San Juan Bay NEP" = "San Juan Bay Estuary Partnership",
"Albemarle-Pamlico Estuary Program" = "Albemarle-Pamlico National Estuary Program",
"Barnegat Bay Estuary Program" = "Barnegat Bay Partnership",
"Delaware Estuary Program" = "Partnership for the Delaware Estuary",
"Delaware Inland Bays Estuary Program" = "Delaware Center for the Inland Bays",
"New York-New Jersey Harbor Estuary Program" = "New York - New Jersey Harbor Estuary Program",
"Massachusetts Bays NEP" = "Massachusetts Bays National Estuary Program",
"Barataria-Terrebonne Estuary Program" = "Barataria-Terrebonne National Estuary Program",
"Morro Bay Estuary Program" = "Morro Bay National Estuary Program")) %>%
group_by(NEP.Name, Year, Activity) %>%
mutate(acres = as.numeric(Acres)) %>%
summarise(totacres = sum(acres))
hab_sum_sf <- left_join(nep_centroid, hab_sum, by = c('NEP_NAME' = 'NEP.Name'))
hab_sum_sf <- spread(hab_sum_sf, Activity, totacres) %>%
mutate(Total = rowSums(cbind(Unspecified,Enhancement,Establishment,Maintenance,
Protection,Reestablishment,Rehabilitation), na.rm=TRUE))
nep_coord <- st_coordinates(hab_sum_sf)
hab_sum_sf <- cbind(hab_sum_sf,nep_coord) %>%
as.data.frame()
hab_sum_sf[is.na(hab_sum_sf)] <- 0
colors <- brewer.pal(7, 'Set1')
hab_sum_sf %>%
leaflet(width = "100%") %>%
addProviderTiles(provider = "CartoDB.Positron") %>%
addMinicharts(hab_sum_sf$X, hab_sum_sf$Y,
type = "pie",
layerId = hab_sum_sf$NEP_NAME,
chartdata = hab_sum_sf[,c(10:16)],
time = hab_sum_sf$Year,
colorPalette = colors,
width = 60 * sqrt(hab_sum_sf$Total) / sqrt(max(hab_sum_sf$Total)), transitionTime = 0)